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ABSTRACT 


A  finite  element  analysis  of  stationary  and  propagating  cracks  in 
the  presence  of  inertia  forces  is  presented.  An  extension  of  the  J-integral 
approach  is  employed.  To  model  a  propagating  crack,  a  conceptually  simple  yet 
effective  technique  has  been  developed.  The  new  crack  propagation  scheme 
eliminates  the  difficulties  associated  with  the  use  of  moving  singular 
elements. 


INTRODUCTION 


It  is  generally  accepted  that  a  crack  arrest  methodology  based  on  a 
dynamic  view  of  crack  propagation  and  arrest  is  more  fundamental  than  the 
quasi-static  approach  to  the  problem  (1).  Specific* lly,  when  inertia  forces, 
stress-wave  reflections,  and  rate-dependent  fractur  processes  are  dominant, 
quasi-static  assumptions  will  generally  underestimi  te  the  true  crack  driving 
force.  In  many  applications,  such  as  nuclear  prest are  vessel  design  and 
structures  subjected  to  Impact  loading,  dynamic  ef '  ects  can  be  important.  In 
these  cases  analytical  models  based  on  the  quasi-static  assumption  may  lead  to 
erroneous  conclusions. 
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Inclusion  of  dynamic  features  in  an  analytical  model  undoubtedly 
leads  to  complications.  In  recent  years  the  finite  element  method  has  emerged 
as  an  important  tool  that  can  be  used  to  resolve  at  least  some  of  the  mathe¬ 
matical  difficulties  associated  with  the  problem.  Considerable  progress  has 
also  been  made  toward  understanding  some  of  the  basic  concepts. 

The  purpose  of  the  paper  is  to  present  the  salient  features  of  a 
recently  developed  finite  element  dynamic  crack  propagation  modeling  tech¬ 
nique.  Results  of  an  elasto-dynamic  analyses  will  be  presented  for  both 
stationary  and  running  cracks.  Comparisons  will  be  made  with  available 
experimental  results.  Through  the  analyses  of  test  specimens,  it  will  be 
demonstrated  how  this  analytical  model  can  be  used  to  acquire  a  better 
understanding  of  the  dynamic  crack  propagation  phenomenon. 

Dynamic  Analysis  and  Numerical  Integration 

The  familiar  finite  element  discretized  version  of  the  equations  of 
motion  are  simply  written  in  the  following  forms 

[M]  {x}  +  [C]  { X}  +  [K]  { X}  -  { R}  (1) 

where 

[M]  -  mass  matrix 

[C]  •  damping  matrix 

[K]  ■  stiffness  matrix 

{ r}  -  the  external  load  vector 

{x} ,  {x},  {x}  ■  the  displacement,  velocity,  and  acceleration  vectors 
respectively. 

There  are  several  methods  that  could  be  used  to  perform  the  direct  numerical 
integrations  of  the  equations.  The  method  selected  for  this  work  was  the 
Newmark  implicit  scheme  [2].  For  this  method,  the  following  approximations 
are  used: 
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•W  ■  K  +  \  +  ixt  +  lc]  4t  <2) 

xc  +  it  "  xt  +  K  4t  +  1(1/2  *  “>  xt  +  “  xt  +  it14t2  <3> 

where  a  and  6  are  parameters  that  can  be  varied  for  accuracy  and  stability 
while  At  is  the  time  step.  From  the  above  equations  it  is  possible  to  write 
the  equations  of  equilibrium  for  time,  t  +  At,  in  terms  of  displacements,  ve¬ 
locities,  and  accelerations  at  time,  t. 

[[K]  +  ajM]  +  djlC]]  {Xt  +  At}  -  {Rt  +At}  +  [M]  (ao{Xt}  +  a^}  +  a3(Xt}) 

+  [CKa^X,.}  +  a4lXt}  +  a5{Xt})  (4) 

where 

*o  •  “2  !  a2  ■  5Jt  !  *3  ■  -  1 

aAt 

a4  "  a  “  1  ;  a5  "  (o  "  2)  5  a6  "  At  (1-<S)  (5) 

a^  -  6At  . 

Solving  Equation  (4)  yields  (x  +  }  whereupon  the  accelerations  and  veloci¬ 

ties  at  time,  t  +  At,  can  be  calculated  using  Equations  (2)  and  (3). 

Elasto-dynamic  Analysis  of  Stationary  Cracks 

Consider  the  problem  of  determining  the  stress  intensity  factor  for 
a  stationary  crack  in  a  structure  subjected  to  dynamic  loadings.  Several  in¬ 
vestigators  have  solved  this  problem  by  employing  singular  elements  around  the 
crack  tip  [3-5].  While  this  approach  has  proved  to  be  successful,  it  will 
later  require  special  considerations  when  the  crack  is  propagating.  An  alter¬ 
native  to  the  singular  element  approach  is  to  derive  the  stress  intensity  fac¬ 
tors  from  a  path  independent  J-integral  [6].  After  accounting  for  the  inertia 
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effects,  the  rate  of  energy  release  per  unit  of  crack  advance  In  the  direction 
of  the  crack  Xj  is  defined  as: 


J  '  [+r  lMdX2  '  °lJnJ  3^  dsl  +  “  •  <6) 

s 

where  the  first  integral  is  the  conventional  J-integral  over  an  arbitrary  path 
surrounding  the  crack  tip  and  the  second  integral  is  an  integral  over  the 
area,  A,  enclosed  by  the  path,  T  +  rg;  see  Figure  1. 

An  application  of  this  J-integral  is  shown  in  the  following  example 
in  which  a  centrally  cracked  plate  in  plane  strain  is  impulsively  loaded  h'  a 
uniform  stress;  see  Figures  2  and  3.  The  finite  element  model  employed  309 
nodes  and  90  eight  noded  quadratic  isoparametric  elements;  see  Figure  4.  The 
material  is  linear  elastic  with  a  Young’s  modulus  of  200  GPa,  Poisson's  ratio 
of  0.3  and  density  of  5  g/cm^.  Newmark’s  implicit  time  integration  scheme  was 
used  (a  ■  0.25  and  6  -  0.5)  with  a  time  step  of  0.15  microseconds.  A  consis¬ 
tent  mass  formulation  was  used  for  this  analysis. 

The  results  of  the  present  analysis  are  shown  in  Figure  5  along  with 
those  of  other  investigators  for  comparison  [7,8].  As  the  figure  shows,  the 
present  analysis  is  in  excellent  agreement  with  the  other  results.  It  should 
be  mentioned  that  the  path  independence  of  this  J-integral  has  been  previously 
demonstrated  [7], 


An  Analysis  of  Running  Cracks 

A  technique  has  been  developed  that  allows  for  the  analysis  of  run¬ 
ning  cracks  without  the  need  for  mesh  adjustments  or  iterations.  This  tech¬ 
nique  is  based  on  earlier  crack  propagation  investigations  using  a  finite  dif¬ 
ference  based  analysis  [9].  It  begins  by  subdividing  the  element  immediately 
ahead  of  the  crack  tip  into  what  can  be  thought  of  as  subelements,  as  shown  in 
Figure  6,  (in  this  case,  for  a  mesh  composed  of  four  noded  linear  elements). 
During  propagation,  the  crack  tip  will  be,  in  theory,  allowed  to  move  in 
discrete  jumps  along  the  crack  plane  through  these  subelements;  e.g.,  the 
crack  tip  will  move  from  Point  "1"  to  Point  "2"  in  one  jump  and  later  move 


FIGURE  1.  COORDINATE  SYSTEM  AND  CURVES  r  and  T, 

cr(t) 


♦  ♦♦♦♦♦ 


a(t) 


FIGURE  2.  GEOMETRY  FOR  STATIONARY  CRACK  ANALYSIS 
a  ■  0.24  cm. 


FIGURE  4.  FINITE  ELEMENT  MESH  AND  J-INTEGRAL  PATH  FOR 
A  STATIONARY  CRACK  PROBLEM 


Normalized  Dynamic  Stress  Intensity  Foe  tor, 

Kj /(<ryira  ) 
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FIGURE  5.  DYNAMIC  STRESS  INTENSITY  FACTOR  VERSUS  TIME  FOR 
AN  IMPULSIVELY  LOADED  CENTER-CRACKED  PANEL  WITH 
A  STATIONARY  CRACK 
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from  Point  "2"  to  Point  "3"  and  so  on,  as  the  analysis  dictates.  The  stress 
intensity  factor  at  any  crack  location  will  be  determined  by  a  J-integral 
evaluation,  Equation  (6). 

Consider  the  situation  with  the  crack  tip  at  Point  "l".  The  crack 
velocity,  V,  will  be  estimated  by  counting  the  number  of  time  steps,  n,  which 
are  required  until  the  stress  intensity  factor,  K,  is  equal  to  the  fracture 
toughness  and  performing  the  simple  calculation: 


where  At  is  the  time  step  size.  This  algorithm  is  repeated  as  the  crack  tip 
moves  from  Point  "l"  to  "2",  "2"  to  "3",  and  so  on. 

This  algorithm  has  the  advantage  of  allowing  the  calculation  of 
pseudo  crack  velocities  after  each  time  step.  This  pseudo  velocity  is  equal 
to  the  actual  crack  velocity  when  the  fracture  criterion 


K  -  Kd  (V,  T) 

is  satisfied,  where  the  fracture  toughness  can  be  a  function  of  the  crack 
velocity,  V,  and  material  temperature,  T.  Since  the  crack  velocity  is  calcu¬ 
lated  first  and  then  the  fracture  toughness,  there  are  no  problems  with  using 
a  velocity  independent  fracture  toughness  relation.  This  type  of  criterion 
would  cause  problems  if  the  crack  was  propagating,  i.e.,  K  ■  KD,  and  the 
inverse  of  the  KD  relation  was  used  to  determine  the  crack  velocity  as  is  done 
in  some  codes. 

It  was  then  necessary  to  decide  the  number  of  subdivisions  to  make 
in  each  element.  This  was  accomplished  by  determining  the  maximum  crack 
velocity  to  be  allowed  in  the  analysis,  Vmax.  In  this  study,  Vmax  was  taken 
as  the  bar  wave  speed,  /E/p  .  Once  Vmax  is  set,  the  maximum  distance  the  crack 
can  propagate  per  time  step,  Ax',  is 


Ax'  «  V  •  At 
max 


(8) 
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The  number  of  Ax'  units  in  an  element  length.  Ax,  is,  To  make  the  number 

of  subdivisions  an  integer  value  and  to  prevent  the  crack  from  traveling  more 
than  one  element  length  per  time  step,  the  number  of  subdivisions,  N,  is  taken 
to  be 


N 


Ax 

Ax' 


+  1 


(9) 


where  N  is  a  truncated  integer  value.  Should  the  crack  want  to  propagate  at 
one  subdivision  per  time  step,  its  propagation  speed  will  be  somewhat  less 
than  the  maximum  velocity,  Vmax,  set  previously  due  to  the  truncation.  The 
amount  by  which  the  actual  maximum  crack  velocity  differs  from  Vmax  will 
depend  on  the  value  of  11. 

The  last  detail  was  to  determine  how  the  crack  tip  would  be  placed 
at  the  subdivision  lines.  This  was  conceptually  performed  by  placing  a  force 
on  the  element  to  which  the  crack  tip  is  adjacent.  For  a  four-noded  element, 
the  force  is  placed  on  the  one  node  on  the  crack  plane  behind  the  crack  tip 
and  for  an  eight-noded  element  nodal,  forces  would  be  placed  on  the  two  nodes 
as  shown  in  Figure  7.  The  forces  were  postulated  to  be  linearly  related  to 
the  crack  tip  location  by  the  following  equation: 


where  F^  is  the  force  at  node  "i”,  F0^  is  the  nodal  force  at  node  "i"  just 
prior  to  the  node  release,  as  shown  in  Figure  7,  "a"  is  the  crack  length  in 
the  element  which  the  crack  tip  is  in,  and  "Ax"  is  the  length  of  that 
element.  In  the  present  study,  eight-noded  quadratic  isoparametric  elements 
were  used.  The  midside  node  was  released  simultaneously  with  the  trailing 
vertex  node.  The  force  history  of  each  node  followed  the  prescription  given 
in  Equation  (10).  Both  nodes  were  released  simultaneously  to  avoid  possible 
problems  discovered  in  another  investigation  [10]. 

A  summary  of  the  above  algorithm  for  a  quasi-statically  initiated 
event  is  given  in  Figure  8. 

In  the  algorithm  just  described,  the  location  of  the  "crack  tip"  is 
obviously  not  actually  known  when  there  are  forces  on  the  crack  face.  The 


TT-. g  ■  g 


Crock  Tip 


■Crock  Tip 


FIGURE  7.  PLACEMENT  OF  NODAL  FORCES  FOR  (a)  A  FOUR-NODED 
ELEMENT  AND  (b)  AN  EIGHT-N0D2D  ELEMENT. 


YES 


FIGURE  8.  FLOW  CHART  OF  PROPAGATION  ALGORITHM  FOR  A 
QUASI-3TATICALLY  INITIATED  DYNAMIC  EVENT 
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only  time  when  the  location  of  the  crack  tip  is  unambiguous  is  when  there  are 
no  forces  on  the  crack  face.  But  in  order  to  avoid  regeneration  of  the  mesh 
(which  would  be  necessary  to  place  a  node  at  the  desired  crack  tip  location) 
the  interpretation  used  in  the  above  algorithm  was  adopted.  A  similar  inter¬ 
pretation  has  also  been  used  by  others  [10,11,12]. 

The  applicability  of  this  interpretation  can  be  somewhat  tested  by 
performing  a  constant  velocity  crack  analysis.  The  problem  chosen  was  a  cen¬ 
trally  cracked  square  plate  of  40  mm  x  40  mm  with  an  initial  crack  length  of 
0.2  a/w,  where  w  is  the  panel  width.  The  panel  was  initially  loaded  with  a 
uniform  stress  in  the  direction  perpendicular  to  the  crack.  The  properties  of 
the  plate  were  E  ■  7716  kgf/mm^;  v  -  0,286;  and  p  =■  2.5  x  10”*^  kgf  *  s^/mm^. 
The  crack  was  propagated  at  a  velocity  of  0.2  of  the  shear  wave  speed  of  the 
material  which  was  Cg  «  3.461  x  10^  mm/sec.  This  problem  is  similar  to  that 
addressed  by  Broberg  [13],  except  that  Broberg  treated  the  crack  as  opening 
from  a  zero  initial  length.  The  mesh  used  for  the  problem  is  shown  in  Figure 
9.  The  model  employed  213  nodes  and  60  eight-noded  isoparametric  elements. 

The  time  integration  employed  a  time  step  of  0.2887  microsecond  with  a  »  0.25 
and  6  -0.5  The  results  of  the  analysis  is  shown  in  Figure  10  along  with  the 
Broberg  solution  and  solutions  of  Atluri  [14].  The  figure  shows  that  the  two 
computed  solutions  are  very  similar  and  they  converge  to  the  Broberg  solution 
after  the  initial  transient  conditions  have  past.  Hence,  very  reasonable 
estimates  of  the  stress  intensity  factor  with  time  can  be  obtained  for  the 
problem  of  a  crack  with  prescribed  velocity. 

The  inverse  of  the  problem  performed  above  is  perhaps  of  even  more 
importance;  i.e.,  a  problem  in  which  the  dynamic  fracture  toughness  is  speci¬ 
fied  and  the  crack  length  time  history  is  sought.  The  ability  of  the  proposed 
procedure  to  perform  such  an  analysis  was  tested  by  comparing  computed  results 
with  experimental  data  for  a  quasi-statically  initiated  crack  in  a  4340  steel 
three-point  bend  (dynamic  tear)  specimen;  see  Figure  11.  The  finite  element 
mesh  used  for  the  analysis  used  314  nodes  and  91  eight-noded  isoparametric 
elements;  see  Figure  12.  A  previously  determined  dynamic  fracture  toughness 
relation  given  by  Kjq  ■  65  +  0.44  V  was  used  [15]. 

The  results  of  the  analysis  given  in  Figure  13  show  excellent  agree¬ 
ment  with  experimental  results. 
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FIGURE  10.  STRESS  INTENSITY  FACTOR  FOR  A  CRACK  STARTING  AT  aG/W  -  0.2 
AND  PROPAGATING  AT  A  CONSTANT  VELOCITY  OF  V/C  ■  0.2. 
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FIGURE  13.  A  COMPARISON  OF  ANALYSIS  AND  EXPERIMENT  FOR  A 

QUASI-STATICALLY  INITIATED  CRACK  IN  A  4340  DYNAMIC 
TEAR  SPECIMEN. 
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CONCLUSIONS 


A  almple  yet  effective  technique  hea  been  developed  to  perforin 
dynamic  creek  propagation  problama.  The  technique  does  not  require  uaing 
aingular  elernta  or  updating  the  finite  element  taeah  during  crack  propaga- 
tlon.  Alao,  the  dynamic  fracture  criterion  need  not  be  velocity  dependent. 
Good  agreement  haa  been  obtained  between  analytical  rraulta  and  experimental 
data. 
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